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Abstract 

We study the formation of vortices in a dilute Bose-Einstein condensate confined in a rotating 
anisotropic trap. We find that the number of vortices and angular momentum attained by the 
condensate depends upon the rotation history of the trap and on the number of vortices present in 
the condensate initially. A simplified model based on hydrodynamic equations is developed, and 
used to explain this effect in terms of a shift in the resonance frequency of the quadrupole mode of 
the condensate in the presence of a vortex lattice. Differences between the spin-up and spin-down 
response of the condensate are found, demonstrating hysteresis phenomena in this system. 

PACS numbers: 03.75.Lm, 03.75.Kk, 67.40.Vs 



I. INTRODUCTION 



The realization of Bose-Einstein condensation (BEC) in ultracold gases has provided a 
powerful new system in the study of superfiuids, combining better experimental control 
and theoretical tractability than liquid helium. A particularly striking example concerns 
the response of the fluid to rotation. Due to the irrotationality property of superfiuids, 
application of a rotating perturbation results in the creation of a lattice of quantized vortex 
lines, a process that was originally observed in He-II Q| . Analogous experiments in ultracold 
gases have found the formation of similar structures in Bose condensates 
and recently in a fermionic gas that is superfluid due to Cooper pairing between atoms [9|. 

An interesting question raised by these experiments concerns the angular velocity of the 
rotating potential required to nucleate vortices. Since the energy per particle of a condensate 
containing a vortex, E^, exceeds that of the ground state, Eq, by transforming to a frame 
rotating with angular velocity fl it is straightforward to show that vortices are energetically 
favored if > (E^ — En)/h [lO|. However, higher rotation rates are generally required to 
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nucleate vortices in experiments [^j S ly, |8|] , which reflects the existence of an energy barrier 



to vortex formation 



This can be overcome by excitation of collective modes localized 



at the surface of the condensate [1^ 

The presence of two different critical angular velocities implies the existence of metastable 
states, where vortices are energetically favorable but dynamically are unable to form, and 



in turn the possibility of hysteresis phenomena |12l . Il^ . Hysteresis has long been observed 
in rotating superfluid helium. For example, experiments have found a significant difference 
between the number of vortices obtained when the fluid is spun up to when it is spun 
down Q|. Reproducible hysteresis loops were measured by minimizing the mechanical and 
thermal noise of the apparatus Q|. Jones et al. [17J attempted to explain this effect 
by using a principle of local momentum conservation. Fortunately atomic Bose-Einstein 
condensates are less noisily coupled to their environment than rotating liquid helium, and, 
unlike liquid helium, can be studied using a model that has great predictive power, namely 
the Gross-Pitaevskii (GP) equation. 

In this article we discuss a particular example of hysteresis, where the amount of cir- 
culation already present in a Bose condensate influences the final angular momentum and 
number of vortices that can be attained. This effect has been exploited in the experiment of 
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Ref. |l8j to rotate the condensate to high angular velocities. The aim of this 
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plore this issue theoretically in more detail. Following previous studies 
we solve the GP equation for the condensate wavefunction \E'(r,i(:) to study the vortex for- 
mation process in a rotating anisotropic trap. However, we go further by considering cases 
where the trap rotation frequency is changed during the simulation, as well as when vor- 
tices are already present in the condensate. Moreover, new insight is gained by comparing 
our findings to the results of a hydrodynamic model 25] of the condensate, which yields a 
resonance in the angular momentum transfer associated with excitation of the quadrupole 
collective mode. The presence of a vortex lattice then has the effect of shifting this mode 
frequency and therefore the resonance to higher values. As a further example of hysteresis 
we consider the case where the trap rotation frequency is slowly ramped up and then down, 
demonstrating a difference in the trajectories of the angular momentum between the two 
cases. 

We note that the discussion in this paper is restricted to very low temperatures, where 
the thermal cloud is insignificant and the condensate dynamics can be accurately modelled 
with the GP equation. An interesting issue concerns the effect of finite temperatures on 
the phenomena described here. However, to address this question requires one to treat the 
dynamics of both the condensate and thermal cloud consistently, since the thermal cloud also 
responds to a rotating potential by spinning up during a timescale related to the frequency 
of collisions between the atoms [26]. Treatment of the coupled dynamics is conceptually rich 
and numerically intensive, and is outside the scope of this work. 

II. VORTEX FORMATION 

The GP equation is solved numerically for a harmonically trapped condensate in 2D, 
which corresponds to the case where the axial trap frequency is much larger than the mean 
field interaction energy, such that motion along the axial direction is frozen. Similar 2D 
studies [2^ 2^ have previously found good quahtative, and in some cases quantitative, 
agreement with experimental data even when this condition is not satisfied [3,0], confirming 
that much of the crucial physics of the vortex formation process is captured in 2D. 

Angular momentum is imparted to the condensate by rotating an elliptically-deformed 
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harmonic trap, which is represented by the potential 

\/(r) = ^mc.l[(l + ey2 + (l-ey2], (1) 

where a rotational transformation of the {x, y, z) Cartesian coordinate system is used, such 
that x' = xcos{Qt) +ysm{Qt) and y' = —xsin{Qt) + ycos{Qt) correspond to rotation of the 
trap at frequency Q. Dimensionless units are also used, with the units of length, time and 
energy given by {h/{muj±)Y^'^, and hw± respectively. The GP equation then becomes 



4^ = [-V^ + (1 + e)x" + (1 - e)y''] + g^'^ ^. 



(2) 



Mean field interactions are represented hj g = A-KN'a, where A^' is the number of atoms per 
unit length in the axial z direction, while a denotes the s-wave scattering length. Throughout 
this paper we will use g = 450, although the results are expected to be generally apphcable 
to other interaction strengths. 

To provide the initial condition for the simulation, Eq. Q is numerically propagated 
in imaginary time with e = and = 0, such that the wavefunction converges to the 
condensate ground state without vortices. In order to model vortex formation the simulation 
is then run in real time with a rotating elliptical trap (e = 0.1, > 0), where the deformation 
e is switched on abruptly at t = 0, rather than being turned on gradually. 

Fig. [T] illustrates the subsequent time-dependent response of the condensate by plotting 
the mean angular momentum (L^) = j dr with = i[ydx — xdy]. HQ = 0.78 (curve 

a) one sees that the angular momentum increases initially, undergoing large oscillations 
before settling to an almost constant value for t > 250, similar to behavior found in previous 
studies as wen a. expense.. J TKe i„,«a. oscUaUon. co„espo„d *o v^iatio^s 

in the quadrupolar deformation of the condensate which subsequently diminish as vortices 
enter the condensate. The vortices then continue to undergo complicated djTiamics, albeit in 
such a way that the angular momentum remains in quasi-equilibrium. However, the vortices 



do not crystal 



dissipation 
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ize to 
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'orm a stable, ordered lattice, which would require the inclusion of 



27L l2q or integration of the GP equation over much longer timescales 



If the trap rotation rate is significantly smaller or larger than Q = 0.78 then the con- 
densate angular momentum displays a periodic time dependence (see Fig. ^ curve b, for 
Q = 0.9). This behavior corresponds to an oscillation of the quadrupolar deformation of the 
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FIG. 1: Angular velocity, (L^) (in units of h) as a function of time, t (units of ujJ_ ) for e = 0.1 
and no vortices present in the condensate initially. The solid lines plot the results when the initial 
trap rotation frequency 0, (in units of uj±) is maintained up to t = 800, with ^1 = 0.78 (curve 
a, black) and = 0.9 (curve b, grey). The broken lines (curves c and d) represent cases where 
17 = 0.7 up to t = 400 (marked by a vertical line), after which the rotation frequency is changed to 
= 0.78 (c, dashed) or = 0.9 (d, dotted). Note that all of the figures in this paper are plotted 
in dimensionless units. 

condensate, 6 = {y'"^ — a;'^)/ (y'^ + x'^). The angular momentum is found to be close to the 
value (Lz) = if the momentum of inertia is given by the expression for a irrotational 
superfluid, 6 = + y'^) j3, 1^- Hence the condensate remains irrotational and no 

vortices are nucleated, in contrast to the case where Q = 0.78. Note that the oscillations 
in the quadrupolar deformation arise due to the sudden imposition of the rotating trap's 
deformation, and their amplitude is lower if e is turned on gradually. 

In the two examples given, the trap rotation rate was maintained at a constant value 
throughout the simulation. If, instead, the rotation rate is changed to another value midway 
through the run, then very different results are found. Curves c and d of Fig. ^ represent 
simulations where for < t < 400 the trap rotates at = 0.7, resulting in vortex nucleation. 
At t = 400 the rotation rate is abruptly changed to either Q = 0.78 (curve c) or Q = 0.9 
(curve d). In the former case the angular momentum increases only slightly, and remains 
much smaller ((L^) ^ 10) than in the case where Q = 0.78 from the beginning ((L^) ^ 18). 
For Q = 0.9 (curve d), however, the angular momentum jumps to a much larger value 
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(corresponding to a major increase in the number of vortices), which is in stark contrast 
with the lack of vortices when the condensate is rotated at this frequency initially (curve b). 
These examples demonstrate that the vortex formation process strongly depends upon the 
"rotation history" of the condensate. 

To explore this question further, we have also conducted simulations where the trap rota- 
tion frequency remains constant throughout each run, but the initial state already contains 
vortices. These states can be obtained by solving the GP equation for a non-inertial frame 
of reference which rotates at constant angular velocity Qq, where t is now imaginary 
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^{-V' + x' + y')+g\^\'-QoL 



^. (3) 



Hereafter we will use Qq to distinguish the rotation of the reference frame (in which we 
solve the GP equation in imaginary time) from the (real time) trap rotation rate Q. For 
sufficiently large Qq, imaginary time propagation leads to the appearance of vortices, which 
eventually settle into an ordered lattice to yield the stationary solution for this rotating 
frame. 

The wavefunction thus found can be used as the initial condition for Eq. Q, which is 
integrated in real time as in the previous section. Fig. |21 shows the resulting time evolution 
of the angular momentum with Q = 0.78 for different initial states. The black solid line 
(curve e) corresponds to the same run represented by curve a in Fig.^ where the condensate 
contains no vortices initially {Qq = 0). The grey line (curve f), by way of contrast, is for 
Qq = 0.5, where the initial state contains six vortices and possesses an angular momentum 
of (Lz) = 3.64. The subsequent rotation of the trap at = 0.78 leads to a final angular 
momentum of (L^) ~ 11, which is only around 60 % of the value attained when no vortices 
are present initially. If higher values of Qq (corresponding to more initial vortices) are used, 
this deficit in the angular momentum becomes even larger. Indeed, if Qq = 0.78, the angular 
momentum remains almost constant throughout the evolution in real time. 

These differences in angular momentum refiect similar variations in the number of vortices 
present in the condensate. This is illustrated in Fig. El where (a)-(c) show snapshots of the 
condensate density, |\E'(r,t)p, at various times for i7o = and Q = 0.78. The corresponding 
case for Qq = 0.6 is shown in Figs. El (d)-(f); vortices are already present initially, but the 
final number of vortices at t = 400 is less than for fio = (comparing Figs. El (c) and (f)). 

To help explain these differences, in Fig. HI we plot the angular momentum as a function 
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FIG. 2: Angular momentum as a function of time for a condensate in a rotating trap with e = 
0.1 and ri = 0.78. The various simulations use different initial conditions, which correspond to 
stationary states for an isotropic condensate in a frame rotating with angular velocity J7o- Plotted 
are = (curve e, solid black line), J7o = 0.5 (curve f, grey), J7o = 0.6 (curve g, dashed) and 
Clo = 0.78 (curve h, dotted). 
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FIG. 3: Condensate density |^'(r,t)p with a trap rotating at angular velocity O = 0.78 for (a) 
t = 100, (b) t = 200, and (c) t = 400, with Oq = 0. Panels (d)-(f) show a similar simulation but 
with rio = 0-6, for (d) t = 50, (e) t = 200, (f) t = 400. The holes in the density profiles correspond 
to vortices. 

of the trap rotation frequency for = and Qq — 0.5. If there are no vortices and (L^) 
has a periodic time dependence then the plotted value is the peak value between t — 
and t — 400. If vortices are present the value of (L^) is taken a,t t — 400, at which 
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point the angular momentum has generally plateaued at an approximately constant value. 
The filled circles represent the results for the initial condition fig = 0, and it is apparent 
that appreciable vortex formation (and hence angular momentum transfer) occurs in the 
approximate range 0.68 < Q < 0.88. Furthermore, within this range the angular momentum 
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rises with increasing f2, in agreement with previous experimental l6|| and theoretical 
studies, where a similar dependence was found when plotting the final number of vortices 
against Q. 




FIG. 4: Maximum angular momentum attained by a condensate in an elliptical trap rotating with 
frequency 0, (e = 0.1). The solid and open circles show the results for f^o = and l^o = 0-5 
respectively. The lines plot the maximum angular momentum found by solving the hydrodynamic 
model (IHira for ^0 = (solid) and Oq = 0.5 (dashed). 



The angular momenta attained at t = 400 when the initial condition is Qq = 0.5 are 
plotted in Fig.Hlwith open circles [3l|. We see that, compared to the f2o = case, appreciable 
increases in (L^) tend to take place at higher Q . Hence, nucleation of additional vortices 
occurs at higher rotation frequencies when vortices are already present in the condensate. 
This effect accounts for the differences in angular momentum apparent in Figs. Q and El 



III. HYDRODYNAMIC MODEL 

The numerical solutions of the Gross-Pitaevskii equation discussed in the previous section 
accurately treat the dynamics of the condensate at very low temperatures. However, it is 
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also instructive to consider a simplified, approximate model that is much less numerically 
intensive to solve than the GP equation. This allows us to more easily explore parameter 
space, as well as providing new insight into the phenomena observed so far. Our simplified 
model is based upon the equations of rotational hydrodynamics, which provide a description 
of the condensate in the Thomas-Fermi (large g) regime 

— + V-(nt;) = 0, (4) 

3v 

— + {vV)v + ViV + gn) = 0. (5) 



Following 



we solve Eqs. (@)) and (EI) to describe the dynamics by employing the following 



ansatz for the density and velocity 

n{r) = ao + a^x^ + ayy"^ + a^yXy, (6) 

(r) = rio X r + Vib^x^ + byv'^ + b^yXy), (7) 

where Oj, bi and fio are time-dependent parameters. Note that the velocity ^ includes 
both rotational and irrotational components, where the former assumes that the combined 
velocity field of the vortices in the lattice approximates a solid body rotation. Substituting 
(jni) and ((Zj) into (jl} and © yields a set of differential equations for each parameter, which 
are integrated in time using a fourth-order Runge-Kutta scheme 

do + 2ao(&^ + by) = 0, (8) 

+ 2axbx + ttxybxy = 0, (9) 

dy - ^ottxy + 2aybx + Qttyby + ttxybxy = 0, (10) 



bx + - (Abl + bly + 1 + + 2gax) = 0, (12) 



dxy - 2QQax + 2Qoay + Aa^yibx + by) + 2(a^. + ay)bxy = 0, (11) 
1 

2' 

by + ^(46J -nl + b% + l- ex + 2gay) = 0, (13) 

bxy + 2(&a- + by)bxy + ey + ga^y = 0, (14) 

do + 2(6, + by)no = 0, (15) 

where = ecos{2Qt) and ey = esin(2r^t), with Q the trap rotation frequency. Once the 
time evolution of these properties are known quantities such as the angular momentum can 
be calculated by integration of (P) and ((Tj). 



The solution of Eqs. (jSlllSp for a rotating elliptical trap yields an angular momentum 
that oscillates in time, similar to the behavior shown in curve b of Fig. The peak angular 
momenta as a function of Q are plotted in Fig. for various values of e and Qq = 0. For 
e = 0.001 the response is symmetrical about a peak near Q = 1/V2. This sharp resonance 
arises since a perturbation rotating with frequency Q is resonant with a surface mode with 
azimuthal quantum number m when u — mQ ~ 0, so that for a quadrupolar trap deformation 
the m = +2 mode should be resonantly excited when fij-es — ^^+2/2 [12]. With increasing e 
the resonance becomes higher and wider as well as more asymmetric, until at e = 0.1 there 
is a steep downward gradient immediately following the peak. Note that the resonance 
peaks are symmetric when Eqs. (jS HT^ are linearized for small departures from equilibrium, 
demonstrating that their asymmetry is a consequence of the nonlinearity of the equations. 
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FIG. 5: Maximum angular momentum (Lz) as a function of trap rotation frequency, fi, solving 
Eqs. (|51IT^ for = 0. The different curves show results for e = 0.001 (solid), e = 0.01 (dashed), 
e = 0.05 (dotted), and e = 0.1 (dot-dashed). The vertical line marks Q = I/V2. 

The results of the hydrodynamic model are compared to those of the GP equation in 
Fig. El The solid line plots the maximum angular momentum as a function of Q for Qq = 0, 
where one sees a close agreement with the GP results. Such correspondence is pleasing, 
but perhaps to be anticipated when vortices are not nucleated, since Eqs. (jHI) and ((7|) are 
expected to be good approximations and the time dependence of (L^) is oscillatory in both 
cases. The agreement is more surprising in the region 0.68 < Q < 0.88, when vortices 
are nucleated and the time dependence departs from oscillatory behavior. It appears that 
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the final, steady angular momentum attained is still close to the peak amplitude of the 
oscillations in the absence of vortices. 

The maximum (L^) for Qq = 0.5 is plotted with the dashed line in Fig.|3J The agreement 
with the GP results is again relatively good, even though the approximation of solid-body 
rotation in ((7j) is not expected to be particularly accurate for the small number of vortices 
present in this case. One sees a shifting of the resonance to higher trap rotation frequencies 
for the larger Qo, which is related to changes in the frequencies of the quadrupole collective 
modes. By linearizing (|HllT3j) for small amplitude oscillations around e = one finds these 
frequencies to be j^] 



This expression provides a simple and useful way to understand the behavior of quadrupole 
mode frequencies with changing rotation. In particular, one sees that the counter- 



while in the presence of a vortex lattice this degeneracy is broken such that uj^2 2 and 
uj-2 for f2o — ^ 1. This general behavior of the m = +2 frequency increasing with Qq is 
reflected in Fig. |31 albeit modified to account for the nonlinearity introduced with e 7^ 0. 

A further consequence of this argument is that to remove angular momentum, and there- 
fore vortices, from the condensate requires the trap to be rotating in the opposite direction 
with an angular velocity approximately half that of the m = —2 mode frequency. We have 
confirmed this by performing simulations with Qq = 0.5 and < 0, and the results are 
presented in Fig. IHl For Q = —0.3 the angular momentum undergoes small oscillations and 
remains close to the original value (L^) ~ 3.64, since this trap rotation rate is sufficiently 
far from the resonance. This is also found to be true when Q = —0.8, albeit with smaller 
amplitude oscillations due to being further from resonance. However, when Q = —0.5 the 
angular momentum rapidly decreases and becomes negative as the vortices that are orig- 
inally present leave the condensate and vortices of opposite sign enter, eventually settling 
into a quasi-equilibrium state. The inset of Fig. El compares the minimum of (Lz) in the 
three simulations to the result of solving the hydrodynamic model (jHllTH|l under the same 
conditions, once again demonstrating reasonable agreement between the two approaches. 
As expected, the peak of the resonance is close to —00^2/'^, where lj_2 — 0.82 from (fTB|l . 
although as before it is skewed towards higher \Q\ due to the large value of e. 




(16) 



propagating m = +2 and m = —2 quadrupole modes are degenerate at uj 



V2 for fio = 0, 
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FIG. 6: Angular momentum as a function of time for = 0.5, and = —0.3 (curve i, solid black 
line), = —0.5 (j, dashed), and 0, = —0.8 (k, grey). Inset: Minimum angular momentum for 
different 0, from GP simulations (points), compared to (solid line) the results of solving Eqs. H8I15|1 
for no = 0.5. 

IV. HYSTERESIS DURING LINEAR RAMPING 

Finally, we again use numerical simulations of the GP equation (j2I) to study the difference 
in the response of the condensate between when it is spun up to when it is spun down. We 
begin with an initial condition containing vortices, such that the wavefunction is a stationary 
solution of Eq. (jH)) with Qq = 0.5. An elliptically-deformed trap potential (e = 0.1) is then 
rotated with a time-dependent angular velocity of the following form 



Hence during the first half of the simulation is linearly ramped up from 0.5 to 0.65, while 
in the second half it is ramped down from 0.65 to 0.5. These correspond to "spin-up" and 
"spin-down" experiments respectively. 

The result is shown in Fig. [71 where the angular momentum is plotted against both Q 
(lower abscissa) and time (upper abscissa). In turn, on the upper abscissa times are shown for 
both the spin-up (black text) and spin-down (grey text) processes, where the corresponding 
plotted curves are represented by the same colors. During the spin- up process (L^) does not 
increase initially, but undergoes oscillatory behavior. This corresponds to the left side of 
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FIG. 7: Hysteresis in the angular momentum, (-Lz), of the condensate wavefunction, initiahy a 
stationary solution of the GP equation Q with = 0.5. The angular velocity of a rotating 
elliptical trap (e = 0.1) follows Eq. ()17() . It is linearly ramped from = 0.5 at t = to = 0.65 
at t = 400 (black line), and subsequently for t > 400, is linearly ramped down such that = 0.5 
at t = 800 (grey line). Note that the abscissa is labelled with $7 (bottom) and t (top), where the 
latter has two different tick labels for spin- up (black) and spin-down (grey). 

Fig.|3J where the rotation frequency is far from the m = +2 quadrupolar resonance. Then, as 
Q approaches resonance angular momentum is transferred to the condensate, until (L^) ~ 6 
at t = 400. During the subsequent spin-down process, however, the angular momentum does 
not follow the path of the spin-up, but remains almost constant with only a slight decrease 
at the end. This reflects the fact that the rotation rate is far from being resonant with the 
m = —2 mode, the excitation of which is required to remove angular momentum from the 
system. Once again this clearly illustrates the importance of the rotation history on the 
attained angular momentum, and hence demonstrates the possibility of observing hysteresis 
phenomena in rotating condensates. 

Finally, we should note that hysteresis behavior was also found by Garcia-Ripoll and 
Perez-Garcia in Ref. Q]. A major difference with respect to the present work is that 
in the time-independent GP equation was solved to find the stationary solutions for 
different rotation rates, while here the dynamics are considered explicitly by solving the time- 
dependent GP equation. Moreover, the results of Ref. Q| are based on the assumption that 
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all collective modes are excited, and vortex formation occurs when the rotation frequency 
Q first becomes resonant with one of the modes. In experiments this would correspond 
to relatively high temperatures where the collective excitations are thermally populated to 
an appreciable extent and the thermal cloud is in rotation (similarly to the experiment of 
Haljan et al. |7]). In contrast, we have been interested in the case where vortex nucleation is 
induced by a rotating elliptical trap, so that the quadrupole mode is predominately excited 
and becomes unstable prior to vortex nucleation. 



V. CONCLUSIONS 



In summary, we have studied the process of vortex formation in a rotating dilute Bose- 
Einstein condensate using numerical simulations of the Gross-Pitaevskii equation. We have 
paid particular attention to the amount of angular momentum transferred to the conden- 
sate for different rotation angular velocities of an elliptically-deformed trapping potential. 
The angular momentum of the condensate (and hence the number of vortices) not only 
depends upon the final rotation rate of the trap, but also upon the history of the rotation. 
Specifically, by initially rotating at one angular velocity, fi, until vortices are nucleated, 
then changing to a second value of fi, one can attain a very different angular momentum 
compared to that achieved when the second VL is retained throughout. Furthermore, we 
have shown that simulations with different initial conditions, corresponding to when vor- 
tices are already present in the condensate, also lead to different final angular momenta, 
similar to behavior found in the experiment of Ref. ^|. Using a hydrodynamic model we 
have demonstrated that this effect is due to a shift of the resonant frequency for the exci- 
tation of quadrupole collective modes when the condensate is already in rotation. We also 
show that the angular momentum response to a linear ramp of the trap rotation frequency 
is radically different depending on whether the ramp is increasing or decreasing, revealing 
the existence of hysteresis phenomena in rotating condensates. 

We conclude by noting that these processes should be readily observable in present exper- 
iments where vortices have been nucleated by stirring the condensate (e.g. Refs. PiSyi), al- 
though the dissipation arising from the presence of a non-condensed cloud may influence the 
angular momentum eventually achieved over longer time scales. Assessing the importance 
of this contribution requires a more sophisticated model that includes flnite temperature 
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effects j35|,|36|, and will be left to future work. 
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